Epidemic Threshold of Susceptible-Infected-Susceptible Model on Complex Networks 
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We demonstrate that the susceptible-infected-susceptible (SIS) model on complex networks can 
have an inactive Griffiths phase characterized by a slow relaxation dynamics. It contrasts with the 
mean field theoretical prediction that the SIS model on complex networks is active at any nonzero 
infection rate. The dynamic fluctuation of infected nodes, ignored in the mean field approach, is 
responsible for the inactive phase. It is proposed that the question whether the epidemic threshold 
of the SIS model on complex networks is zero or not can be resolved by the percolation threshold 
in a model where nodes are occupied in the degree-descending order. Our arguments are supported 
by the numerical studies on scale- free network models. 
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Epidemic spreading is a common phenomenon in net- 
worked systems. Diseases spread from individual to indi- 
vidual through a contact network and computer viruses 
spread through the Internet. Since it has a huge im- 
pact on stability, epidemic spreading on complex net- 
works has been attracting a lot of interests during the 
last decades [l[ . Those studies have focused on both the- 
oretical issues such as noncquilibrium critical phenom- 
ena Q and practical issues such as searching for an ef- 
ficient immunization strategy 0, The SIS model is 
a paradigmatic epidemic spreading model where an in- 
fected individual becomes susceptible (or healthy) at a 
unit rate and infects its susceptible neighbor at a rate A. 
We consider the SIS model on complex networks whose 
degree distribution P(fc) denoting the fraction of nodes 
with degree k is broad [5|. 

Pastor-Satorras and Vespignani proposed a so-called 
heterogeneous mean-field (HMF) theory for complex net- 
works [6|. According to it, the epidemic threshold of 
the SIS model, above which the system is in an active 
phase with a finite density of infected nodes, is given by 
A c = (k)/(k 2 ) with (k n ) = Jdkk n P(k). Specifically, in 
scale-free (SF) networks characterized by P(k) ~ fc~7 
with a degree distribution exponent 7 [5[, A c = for 
7 < 3 while A c > otherwise 0. The HMF theory, 
which becomes exact in the annealed network limit |8|- 
Hoj . turns out to be useful in studying various physical 
problems on complex networks 0] . 

Meanwhile, a recent mean filed (MF) study [ll[ on 
rather realistic quenched networks reports that the epi- 
demic threshold vanishes (A c = 0) in any network with 
diverging maximum degree. It implies that an epidemic 
spreading cannot be prevented on complex networks with 
an unbounded degree distribution. This study attracts 
much interest and is followed by a series of works (l2- 15|. 
The MF theory on quenched networks will be referred to 
as the quenched mean-field (QMF) theory. The discrep- 
ancy between the two MF theories makes it urgent to 
study the SIS model beyond a MF level. 



In this Letter, we challenge the MF approach by tak- 
ing account of the dynamic fluctuation of infected nodes. 
This effect turns out to be crucial in determining whether 
the epidemic threshold of the SIS model on complex net- 
works is zero or not. Wc find that the active phase 
predicted by the QMF theory near A = actually cor- 
responds to the Griffiths phase TB"3] where the den- 
sity of the infected nodes decays to zero more slowly 
than an exponential decay, unless the active nodes in the 
QMF theory form a percolating cluster. It is proposed 
that zero/nonzero epidemic threshold of the SIS model 
is inherited from zero/nonzero percolation threshold in a 
model where nodes are occupied in the degree-descending 
order. Such a specific percolation will be referred to as 
the degree-ordered percolation (DOP). Our argument is 
confirmed in the numerical studies on the (it, u)-flowcr 
model [3] for scale-free networks. We finally apply the 
DOP to survey whether A c of the SIS model should be 
zero or not on random SF networks, which remains un- 
settled in model simulations due to a strong finite-size 
effect 



We begin with a review on the QMF theory for the SIS 
model. Let pi (t) be the infection probability of node i at 
time t. Then the rate equation reads 



dpi (t) 
dt 



= -M*) + (1 -/*(*)) . (1) 



where a%j is an element of the adjacency matrix assigned 
with 1 if there is an edge between nodes i and j or 
otherwise. The first term of the R.H.S. of Eq. ([I) is the 
recovery rate reducing the infection probability and the 
second term is the infection rate given by the product of 
the susceptible probability and the infection trial rate by 
infected neighbors. 

The QMF approach focuses on the linear stability anal- 
ysis of the zero fixed point (p;(0) = for all i) of 
Eq. (JlJ , which corresponds to a configuration of the inac- 
tive phase. It is easy to show that the fixed point becomes 
unstable as soon as AAi > 1 for the largest eigenvalue Ai 



2 



of {ciij}. This leads to the conclusion A^ MF = 1/Ai for 
the epidemic threshold of the QMF theory j2(| • Although 
appealing, it has some controversial points. Most of all, 
it predicts A c = in any network with the diverging 
maximum degree. In an arbitrary graph with the maxi- 
mum degree fc max , the largest eigengenvalue satisfies an 
inequality y^max < Ai < fc max [21] . This gives A^ MF = 
in the fc max — > oo limit. An alternative interpretation of 



\Qmf f ii ows recently in Ref. 14j, which claims that a 
property of the eigenvector corresponding to Ai plays an 
important role in epidemic prevalence. 

As a counterexample to the QMF conclusion, it is in- 
structive to consider a star graph consisting of a hub at 
center and fc max linear chains of length L emanating from 
it. The total number of nodes is N = k mllx L + 1. It is 
straightforward to show that Ai — > k max /y/k max — 1 for 
large L. Hence, A^ MF = in the infinite A: max limit. On 
the other hand, the steady state solution of Eq. (JTJ is 
given by p r oc [2A/(1 + \fl - 4A 2 )] r for 1/Ai < A < 1/2 
and Pr = (2A - 1)/(2A) for A > 1/2, where 1 < r < L 
denotes the distance from the hub. So the epidemic order 
parameter p = lim Lifemax ^ oc £\ Pi/ N is given by 







for A < 1/2 



P = 



A -1/2 for A > 1/2 



(2) 



Namely, A c = 1/2. This example demonstrates that the 
linear stability analysis against the inactive state alone 
is not sufficient in determining the threshold of the SIS 
model on complex networks. 

In order to overcome the shortcoming of the previous 
QMF approach, we suggest to take account of the other 
eigenvectors of the adjacency matrix than the leading one 
with the largest eigenvalue Ai. This approach was also 
taken in Ref. [3]. More importantly, we also suggest 
to take account of the fluctuation of active nodes due to 
the stochasticity in updating each node's state, which is 
ignored in most MF approaches. 

For convenience, we label the nodes in the degree- 
descending order: k\ > fc 2 > • • • > fejv. Recent studies 
show that the n-th largest eigenvalue of the adjacent ma- 
trix of a random network is A„ ~ v^n for large k n (22| 
and the corresponding eigenvector is localized around the 
node n 23J. These findings imply that the steady state 



solution of Eq. ([1) at small values of A consists of local ac- 
tive domains around high degree nodes: Each high degree 
node n behaves like an independent local hub with its 
own activation threshold given by A„ = 1/A„ ~ 
and the size of a local active domain is given by ~ Afc„. 
Independence of the local active domains in the small A 
limit is guaranteed only when higher degree nodes are dis- 
tant enough from each other (this will be clarified later). 
A network with such a property will be referred to as 
the unclustered network. For 1/A„ < A < l/A n+ i, there 
appear the local active domains of size ~ Xki around all 
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FIG. 1. (Color online) Density of infected nodes in the (3, 3)- 
flowers in (a). The largest cluster density in the DOP (see 
text) on the same flowers as in (b). 



nodes i < n. Thus one may expect 



P 



1/A 2 



dk (Xk)P(k) 



(3) 



which yields p ~ A 27 " 3 for a random SF network with 
a degree exponent 7. Note that the refined QMF the- 
ory still predicts that A c = with the order parameter 
exponent ft = 27 — 3. 

A numerical evidence, however, shows Eq. (j3]) is not 
valid. We have performed Monte Carlo simulations for 
the SIS model on SF networks generated by the (u, v)- 
flower model [lj| . It is a deterministic hierarchical model: 
One starts from two nodes connected with an edge (ze- 
roth generation) . Then, every link in a G-th generation is 
replaced with the two u- and w-link-long paths in the next 
(G + l)-th generation. It results in a SF network with 
7 = l+ln(w+u)/ln2 in the G -> 00 limit 0. The (u,v)- 
flower model is particularly useful because one can gener- 
ate an unclustered network easily If u > 1 and v > 1, the 
degree of all nodes is doubled and the distance between 
them becomes farther after each iteration. We used the 
(3, 3)-flower model in simulations which certainly belongs 
to the unclustered network. The numerical data shown 
in Fig.[T](a) strongly suggest a transition at nonzero A c . 
The threshold can be estimated from the peak positions 
of the susceptibility x = N{(p 2 ) — (p) 2 )/(p) where A^ is 
the total number of nodes |l5|. The inset in Fig. Q] shows 
the peak position is extrapolated to A c ~ 0.65(5) in the 
infinite system size limit. 

The origin of the inconsistency is an irreversible fluctu- 
ation within local active domains due to the stochasticity 
in microscopic dynamics. Consider a local domain con- 
sisting of V nodes where each node is in the infected 
state with a probability r and the susceptible state with 
a probability (1 — r). Then, no matter how rare, there 
exists a moment when all nodes recover simultaneously 
by chance. It takes place after a characteristic time 
Ty ~ (1 — r)~ v ~ e rV . Once being recovered, the lo- 
cal domain never returns to an active phase unless cxtcr- 
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FIG. 2. (Color online) p(t) in the (3, 3)-flower of G = 8. 
The data points are the SIS model simulation result. The 
solid curves are obtained from the direct integration of Eq. (|4|) 
with two fitting parameters a± = 0.8 and a-2 = 0.325 for r\k = 
aie a2Xk and with the truncated degree distribution function; 
from bottom to top, P(k) is truncated by k > fcmax/2™ with 
n — 1,2, ..,7. Inset compares the long-time behavior pre- 
dicted in our theory (red solid curve) with the currently avail- 
able numerical data (blue dots). 



generate local domains which behave differently from the 
bulk. Denoting the probability that such a domain of size 
£ is realized as P(£) and the relaxation time therein as 
t(£), a physical quantity / relaxes to its stationary value 
/. as 5f(t) = f(t) - f s ~ /deP(0e-*/ T «>. In the Grif- 
fiths phase, the relaxation dynamics is dominated by rare 
events encoded in the tail of P(£) with long characteristic 
time scales r(£) [24]]. In our case, the slow dynamics is 
originated from the irreversible fluctuation near hubs. 

The previous argument shows that the SIS model on 
the unclustered network can be in the inactive Griffiths 
phase for A^ MF < A < A c with a certain nonzero A c . 
It also provides a hint on the mechanism for the phase 
transition into an active phase. Inside the Griffiths phase, 
local active domains of size ~ Xk are separated [unclus- 
tered network] and metastable. As A increases, the size 
of local domains grow and they begin to overlap each 
other. The active domains become globally stable when 
they will form a percolating giant cluster at a certain 
threshold value of A. That is to say, the epidemic tran- 
sition is triggered by a percolation transition of the local 
active domains. 



nally activated. Considering this effect of the irreversible 
fluctuation, Eq. © should be replaced by 



1/A 2 



dk (\k)P{k)e- t/Txk 



(4) 



Since it decays to zero eventually, the apparent active 
phase implied in Eq. ([3]) is in fact an inactive one. We 
remark that the inactive phase is different from the 
usual one in which the density decays exponentially in 
time. The density in Eq (Ql docs not decay exponen- 
tially fast but extremely slow due to the broad distribu- 
tion of relaxation times. For example, in SF networks 
with P(k) ~ fc -7 , one finds 



p(*)~M 



-(7-2) 



(5) 



Equation ([4]) is numerically tested in the (3, 3)-flowers 
at A = 0.56 < A c , which is shown in Fig. [2] It does 
not decay exponentially in time. Instead, the density 
is trapped in a plateau for a while and then decays to 
another plateau successively. Those plateaus are the ev- 
idence for the existence of the metastable local active 
domains. In the (it, i>)-flower, degrees of nodes are dis- 
cretized as k = 2™ with n = 1, 2, ... So the size Vk ~ Xk of 
the local active domains and their life time T\k ~ e aXk are 
also discrete. This discreteness results in the plateaus. 
As shown, the data in Fig. [2] are fitted to Eq. © very 
well for large t > 10 2 . The dashed curve therein repre- 
sents the overall decay given by Eq. ([5]). 

The slow dynamics given in Eqs. Q and is reminis- 
cent of a relaxation dynamics in the Griffiths phase 0- 
18|,[2J]. In a disordered system, disorder fluctuations may 



Note that the local active domains nucleate around 
high degree nodes in the degree-descending order. So the 
uncovered mechanism leads to the conjecture that the 
unclustered network with A c ^ should have a nonzero 
percolation threshold in the DOP model. Recall that the 
DOP is suggested as a percolation model where nodes are 
occupied in the degree-descending order [25j. A nonzero 
percolation threshold p c for the node occupation proba- 
bility implies that high degree nodes are well separated 
from each other. Hence, in the context of the SIS model, 
one requires a nonzero value of A for the local active do- 
mains (of size ~ Xk) to form a percolating cluster. 

We provide a numerical evidence of our claim. Fig- 
ure [T] (b) shows the the percolation order parameter g, 
the density of nodes in the largest cluster, for the DOP 
in the same (3, 3)-flowers used in (a). The percolation 
threshold p c is clearly nonzero, which is consistent with 
the fact that A c ^ as shown in Fig. [1] (a). As another 
example for the validity of the claim, one may revisit the 
aforementioned /c max -star graph case with large L. In 
this example, one can easily find that p c ^ supporting 
our claim. 

We also consider the opposite case with p c = 0. It is 
achieved only when any finite fraction of occupied nodes 
in the DOP process are connected to each other to form 
a percolation cluster. Those networks with p c = will 
be called the clustered network. In the context of the 
SIS model, the local active domains in the clustered net- 
work form a percolating cluster even in the limit A — > 0. 
Therefore, we expect that the epidemic threshold is zero 
in the clustered network. The steady state density of the 
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FIG. 3. (Color online) Density of infected nodes in the (1, 5)- 
flower in (a). The largest cluster density in the DOP on 
the same flower in (b). The solid lines in (a) has a slope 
21n6/ln2 ~ 5.19 while the dashed line in (b) has a slope 1. 



infected nodes can be expected to scale as 

dkP{k) . (6) 

l/A 2 

In comparison with Eq. ([3]), the factor (Afc) is missing 
because a stable active domain consists mainly of high- 
degree nodes with k > l/A 2 . In scale- free networks with 
P(k) ~ fc -7 , one obtains p ~ A 27-2 . Note that the 
order parameter exponent /3 = 27 — 2 is different from 
/3qmf = 1 obtained from the simple QMF theory where 
only the largest eigenmode was taken into account fl2| . 

A trivial example of the clustered network is the k max - 
star graph with L = 1. The density of the largest cluster 
is given by g = p implying that p c = 0. Independently, 
A c = obviously as discussed in Ref. [HI]. A nontriv- 
ial example is the (u, v) -flowers with u = 1 (or v = 1). 
Recall that if two nodes are connected with an edge in 
a certain generation, then they remain being connected 
afterward when u = 1 (or v = 1). So one can expect that 
high degree nodes are clustered. In Fig. [3][a), we present 
the SIS model simulation results on the (1, 5)-flowers at 
several generations G. As G increases, the data approach 
the theoretical prediction p ~ (A — A c )^ with A c = and 
/3 = 27 - 2 = 2 In 6/ In 2. The DOP property is presented 
in Fig.EJb). The scaling of g ~ p therein indicates p c = 0. 

What is the epidemic threshold in more interesting 
cases such as generic random SF networks? Recently, 
extensive Monte Carlo simulations were performed in 
random SF networks generated from the configuration 
model [IH. However, due to a strong finite-size effect, 
it still remains inconclusive whether A c = or not even 
with simulations of system sizes up to N = 3 x 10 7 . 

Alternatively we investigate the DOP property of the 
configuration model networks. In Fig.|4ja) we present the 
percolation order parameter for the network with 7 = 5. 
As in the (3, 3)-flower, the system undergoes the perco- 
lation transition at a finite threshold. In order to esti- 
mate the percolation threshold precisely, we make use of 



10 



10 



10 



10"' 



10 



1 


1^-— 
















(a) 


1 





-5[ 

0.2 0.4 



0.6 



1 

10 4 


1 


- D 10 5 








10 6 




- A 10 7 




< extra. 






(b) 


4^*1 


1 



3 3.5 4 4.5 5 5.5 6 

y 



FIG. 4. (Color online) (a) The density of the largest DOP 
cluster in random SF networks of the configuration model 
with 7 = 5. The dashed line denotes the percolation tran- 
sition point. (b) Maximum cluster heterogeneity points 
and their extrapolated values. The system sizes are N = 
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a so-called cluster heterogeneity (CH) which denotes the 
number of distinct cluster sizes [27(. It was shown 28 1 
that Pc(N) at which the CH is maximum in networks of 
size N converges to the percolation threshold p c in the 
infinite N limit. In Fig. EJb), we present the numerical 
data for p c (N) at several values of 7 and their extrapo- 
lated values. The percolation threshold is nonzero unless 
a small 7 is considered. Thus, the random SF networks 
therein belong to the unclustered network class. It pro- 
vides an indirect evidence that the epidemic threshold 
could be nonzero on those SF networks. 

In summary, we present a theoretical argument that 
the epidemic threshold of the SIS model on complex net- 
works is nonzero in the unclustered network while it is 
zero in the clustered network. This conclusion is drawn 
by taking account of the effect of the irreversible fluctu- 
ation which was ignored in the QMF theory. The fluc- 
tuation makes a local active domain unstable and leads 
to the Griffiths phase. Numerical simulations performed 
in the deterministic (u, w)-flowers support our argument. 
We suggest that the clustering property of a network can 
be determined by the DOP. By studying the DOP tran- 
sition, the random SF networks are shown to belong to 
the unclustered network unless a small degree exponent 
is considered. It suggests the epidemic threshold in such 
SF networks is nonzero as opposed to the QMF predic- 
tion. Our work raises various interesting questions on the 
critical phenomenon associated with the epidemic tran- 
sition and the DOP transition, which are left for future 
works. 
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